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Abstract 

We have written a Fortran programme BCVEGPY, which is an event generator for the hadronic 
production of the Be meson through the dominant hard subprocess gg — > Bc{B*) + b + c. 
To achieve a compact programme, we have written the amplitude of the subprocess with the 
particle helicity technique and made it as symmetric as possible, by decomposing the gluon self 
couplings and then applying the symmetries. To check the programme, various cross sections 
of the subprocess have been computed numerically and compared with those in the literature. 
BCVEGPY is written in a PYTHIA-compatible format, thus it is easy to implement in PYTHIA. 
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I. INTRODUCTION 



Sc-physics has been attracting an increasing attention recently, due to the experimental 



discovery o: 



16 



the Be meson and theoretical progress 



IG 



11 



12 



171. Since one can collect a high-statistics sample of mesons only at high 
energy hadronic colliders have re-written a generator for the hadronic 

production of B^. The generator is a Fortran program package called BCVEGPY. 

The hadronic production of the Be and B* mesons 12811 has been estimated by' the authors 
of Re..B wi* .a,.en*a«on app^ach and tKe au.ho. of RefsJ B B with 
the so-called 'complete calculation' approach i.e. to compute the production completely at 
the lowest order (a^) in terms of the dominant subprocess of perturbative QCD (pQCD) 
gg —>■ Bc{B*) + c + b. Further note that since mb ^ rric ^ ^qcd, apart from the fact 



the rest of the subprocess is 
These two approaches have 



that c and b quarks combine into BAB*) non-perturbative 
always 'hard' and can be well calculated with pQCDp, |8|, |! 
different advantages but both of them are in the framework of pQCD, and attribute the non- 
perturbative factor in the production to the decay constant, i.e. the non-relativistic wave 
function of the Bf.{B*) at the origin (the former through the fragmentation function and the 
latter directly). Within theoretical uncertainties, the two approaches can agree numerically, 
especially when the component of gluon fragmentation is involved [lo|. The fragmentation 
approach is comparatively simple and can reach to leading logarithm order (LLO) of pQCD, 
but is satisfied only in the case if one is only interested in the produced Bc{B*) itself. 
The complete calculation approach has the great advantage that it retains the information 
about the c and b quark (jets) associated with the B^ meson in the production. From the 
experimental point of view this is a more relevant case. Therefore, we have written the 
hadronic production programme for B^ mesons based on the complete calculation approach 
(full pQCD complete calculation at the lowest order a^). 

Since the non-perturbative factor in the subprocess gg Bc{B*) + ■ ■ ■, i.e. the decay 
constant of the Bc{B*) meson, can be calculated by means of the potential model for heavy 
quark-antiquark systems[ll|, the estimates of the production are full theoretical predictions 
of pQCD without additional experimental input. In view of future experimental studies 
of the Bc{B*) after the CDF discovery (RUN-1, at Tevatron), i.e., concerning the needs of 
experimental feasibility studies for various topics of Bc{B*) meson in hadronic collisions at 
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Tevatron and LHCj29|, we would like to write a paper on the Fortran package BCVEGPY 
with detailed explanation. In addition, we emphasize here that, as explained later, the 
present generator has been formulated in a very different way from those in Refs.p, O], 
and through carefully comparing the results obtained by BCVEGPY with those obtained in 
Refs.0, 01, a solid and independent check of BCVEGPY is well made. 

At LHC the beam luminosity and the production cross section (at such a high energy 
= 14TeV) are so high that the rate of producing Bc{B*) events can be 10*^"^ per year. 
At Tevatron, although the luminosity and production cross section are lower than at LHC, 
the rate of producing Bc{B*) events still is about 10"^"^ per year, only 3 ~ 4 orders of 
magnitude lower than at LHC. The B^ meson has sizable and abundant weak decay modes 
relating directly to c or/and b flavor respectively, so that with so high production such as 
at LHC, even the detecting efficiencies being taken into account, one can reach a statistical 



accuracy of 10 ^ for most of the decay modes 
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1J|. Thus there is an acute need for 



a Bc{B*) event generator in order to be able to perform feasibility studies, and this is why 
we have written the present paper on BCVEGPY. A particularly interesting topic, which 
is worth noting here, is the study of Bg — Bg mixing and CP violation in Bs meson decays 
through B,{B*) ^ The Be meson has a very large branching ratio to decay 

to a Bg meson {Br {Be ^ Bs ■ ■ ■) — several tens percentsjl^. Iisl 1^). and the Bg mesons 



obtained by such B^ decays are tagged precisely at the B^ decay vertex, if the charge of 
the Be meson and/or its decay products are measured. A large-statistics sample of tagged 
Bs mesons at LHC, and even at Tevatron, offers a great potential for the interesting Bs 
physics 0]. 

As pointed out in Ref.jTj, according to pQCD there is another production mechanism by 
a quark pair annihilation subprocess qq Bc{B*) + c + 6. Nevertheless, the 'luminosity' 
of gluons is much higher than that of quarks in pp collisions (LHC) and in pp collisions 
(Tevatron), and there is a suppression factor due to the virtual gluon propagator in the 
annihilation. Therefore the contribution from this mechanism is negligible compared to the 
dominant one. TKe calculations in Rcfs.0 3 3 Q neglected the contnbut.on .on, tl. 
quark pair mechanism, and the BCVEGPY package follows the same approximation. 

In order to make the programme very compact, we write BCVEGPY by applying the 
'helicity technique' to the amplitude of the subprocess. The technique may be traced back 
to the work in Ref . [iS^ . The helicity technique has been developed by the CALKUL 
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collaboration 



la |20| for massless gauge theories. Further development for the massive 



fermion case with Abelian gauge field(s) as well as for the massless fermion case with non- 
AbeUa„ gauge «e.d(s) was doe by seve.a. groupsQ Q H Q. Whe. the KeliCiesfl 

of all the external massless particles are fixed, CALKUL calculates the 'probability rate' in 
the following way. First, Feynman diagrams with precise helicities for the external particles 
are computed one by one for the concerned process with helicity techniques, and a complex 
number for each Feynman diagram is obtained. Then all the obtained complex numbers are 
summed up. Finally the squared modulus of the summed result is taken, averaging over the 
helicities of the particles in initial state and summing over the helicities of the particles in 
final state if an unpolarized case is considered. If a polarized process is considered, one needs 
only to stop the calculation before averaging over the helicities of the particles in initiate 
state and summing over the helicities of the particles in final state. 

The massless spinor technique may, in fact, be applied to the case where massive 
fermion(s) and non-Abelian gauge field(s) are involved in the concerned process, if a suit- 
able 'generalization' could be done together with some further rearrangements. Our present 
subprocess contains non-Abelian gluons and massive fermions. Thus, in order to make use 
of the massless spinor technique for obtaining a compact result, we need a suitable general- 
ization with appropriate rearrangements. We will describe the procedure here below. Our 
strategy for the generalization is to convert the problem into an equivalent 'massless' one 
and to extend the 'symmetries' as much as possible. Then we try to apply the symmetries 
of the converted amplitude and the helicity technique to the present problem so as to make 
the programme compact. According to pQCD, at order there are 36 Feynman diagrams 
for the 'hard' subprocess gg Bc + c + b. To extend the symmetries for the amplitude corre- 
sponding to the 36 diagrams, we neither consider the color factors nor distinguish the flavors 
of the fermion lines at the moment. Then, these diagrams may be grouped into only a few 
typical ones according to the fermion lines and the structures of the contained 7-matrices on 
the lines in the Feynman diagrams, because of the Feynman diagram symmetries. To apply 
the symmetries in writing up the program, we first focus on the numerator of the ampli- 
tude related to each typical fermion line, and deal with the 7-matrices precisely. We then, 
having suitable four-momentum set for the numerator and denominator, respectively imple- 
ment proper numerator factors for the fermion lines: color factors, suitable denominator and 
spinors (corresponding to the external lines) etc. Finally we obtain an exact and full typical 
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fermion line, which appears in Feynman diagrams. When all kinds of typical fermion line 
factors, factors for external lines of gluons and gluon propagators are 'assembled', then the 
full term, corresponding to the Feynman diagram of the amplitude, is achieved. The result- 
ing program is indeed very compact and potentially reduces the execution time significantly. 
In general, to write an amplitude according to helicities, each massive fermion line should 



be decomposed into two light-like spinor lines, but this is not unique |19l. 1201 l21l. l22l l23l l24 1 . 
Here we use the identity f, = \k+){k + | + \k—){k — \ (fc is a light-like four-momentum) to 
simplify the fermion (quark) lines by spinor products. 

To verify that the program BCVEGPY is correct, we check it by taking the same pa- 
rameters as in Ref. ^, compute the cross-section of the subprocess gg Bc{B*) + c + 6 by 
integrating out the unobserved variables numerically, and compare the obtained results with 
those in Ref. 0] carefully. Since BCVEGPY is based on helicity techniques which is totally 
different from those in Ref. ^, the comparison is a very good check both for BCVEGPY and 
for that in Ref. 3. 

Since the mass difference between -Bc[l/ •S'o] and B*[l,^ Si] is small, the B* decays through 
an electromagnetic transition (Ml) B* ^ Be + '~f with an almost 100% branching ratio and 
the photon in the decay is quite soft, so the production of B* can be considered as an 
additional source of Be production (with an additional soft photon). Furthermore, since 

= 0~ for Be, the polarizations of b and/or c-jets will be lost during their hadronization, 
so the polarization effects including those of b and/or c-jets in the production of Be are 
therefore not interesting, so BCVEGPY program calculates the process for the un-polarized 
case only, although theoretically it can work out certain polarization effects due to B* and/or 
b and/or c-jets. 



We write the BCVEGPY package following the format of PYTHIA[25| so that the gener- 
ator could be easily adapted into the PYTHIA environment. In this way our generator can 
be used for generating complete events. To increase the Monte Carlo simulation efficiency 
for high dimensional phase space integration, we set up a switch in BCVEGPY to choose 
if the VEGAS program for obtaining the sampling importance function is used or not. We 
also use several parameters, such as for the maximum differential cross sections etc, to meet 
the needs for the initializations of PYTHIA. 

This paper is organized as follows. In Section II we show how to extend the symmetries 
by: focussing the 7-matrix strings of fermion lines in Feynman diagrams, disregarding color- 
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and numeral factors etc; decomposing the three- and four-gluon couphng vertices; grouping 
the decomposed diagrams; estabhshing the typical one in each group; applying the sym- 
metries and helicity techniques to the problem. In Section III we outline the structure of 
BCVEGPY, explain how to use the programme and test (check) the programme as stated in 
the Introduction. Section IV summarizes the conclusions and future prospects. The details 
about the hehcity functions for the amplitude, polarization vector for B*fSi], routines and 
functions for the helicity amplitude are described in the Appendices. 



II. THE HARD SUBPROCESS 



Based on the factorization of perturbative QCD (pQCD), the hard subprocess play the 
key role of the Be generator. In hadron collisions at high energies, the subprocess gluon- 
gluon fusion gg Bc{B*) + 6 + c is the dominant one in hadronic production of Bc{B*) 
mesons j^. In the following sub-sections, we will show how to deal with this subprocess. 



A. The amplitude for gg Bc{B*) + b + c 

At the lowest order a^, there are totally 36 Feynman diagrams as shown in Figs. 1-5 
for the gluon-gluon fusion process gg —>■ B^ + b + c, so accordingly there are 36 terms in 
the production amplitude. As stated in the Introduction, we group the Feynman diagrams 
into five sets according to the structure of the fermion lines. Here there are five subsets 
as indicated in Figs. 1-5. To write the amplitude corresponding to the Feynman diagrams 
explicitly, let us first focus on writing the color structures of the diagrams separately. In fact 
there are five independent color factors, Cuj, C2ij, C^ij, C^ij and C^ij, where i,j = 1,2,3 
are the color indices of the quarks c and b respectively. They are 

- 1 

N"^ - 1 

C,^, = -\6.,Tr[T'^T'] . (1) 
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FIG. 1: Feynman diagrams that can be directly grouped into the cc subset. Here i and j are the 
color indices of c and b respectively. 

To make the final amplitude compact, we introduce two extra color factors, C'^^j and Cgj^, 
which satisfy 



(2) 



All the color factors in the subprocess can be expressed by the above five independent 
color factors, i.e. in the amplitude all the color factors may be written in terms of these 
five explicitly. To obtain the desired result, the general commutation relation for the color 
factors 

[Ta,Tfc] = ifabcTc (3) 

has been applied to three- or four-gluon vertice of the Feynman diagrams, where fahc is the 
antisymmetric SU{?>) structure constant. 

The terms of the amplitude corresponding to the grouped Feynman diagrams are written 
as below. 

The first group, FiglH 

J (27r)4 [ (g^i + qf,2y {h + ^2 - qc2)^ - mj 
J {2ttP [ {qbi + qb2r [ki + ^2 - qc2r - 



M- 



le 



M- 



1/ 



Ml 



19 



M- 



Ih 



'1 - 4c2 + m 



{h - qc2y - ml ^ 



Mic = CiijUs{qbi)i j 



(2vr)4 



2^2 



Mid = C3ijUsiqbi)i 



(2vr) 



Vs'{qc2) , 

Xp(g) /Ai - + me 

^ {qbi + qb2Y ^ (gel - - m^ ^ 

Vs'{qc2) , 



75 



rf2 



(gw + g62)^ ^ (gel - ^^2)^ - ml 

Vs'{qc2) , 



75- 



Xp(g) il-l^l + ^e 

^'(gw + g62)2^' (gel - fci)2 - mf 2 ' 



f 1 - ge2 + mc . 

{h-qc2f-mf' 
^6^i^^(gw)^ / 

^el - 1^1 - ^2 + m^ 

(gel - fci - - ml 

c^5.-.(g.i)^/^{7.^ 

^el -h-h + ^c 

(gel - A^i - /i;2)^ - ml 

{C2ij - Cuj)Us{qbi)i I l757-^^7^^^75: 



75 1 ^^s'(ge2) 

xp(g) 



,A2 _gcl-|2 + mc Ai 



hi + g62)^ ^ (gel - 

15 \ Vs'{qc2) , 



2 ''I 



^1 + ^2 - ^c2 + mc 



(27r)4 i '"(9,1 + ^,2)2 "'(fci + A;2-ge2)^ 



[Kl + K2) ^ ^ J 



fs'(gc2) 



(c^6. -c^5.-K(g.i)^y (^|7.^ 



xp(g) 



+ mc 



Ibi + qb2Y {ki + k2Y {ki + k2- qdf 



75 



[{ki - k2)ag,iiy + (2/C2 + ki)f,g^a + (-2/ci - k2)^g^a) |fs'(gc2) 



The second group, FigI21 



M2a = (^6.,)^.(g.l)^y (^|75^ 



2b 



M. 



2c 



d^q j h + h- ib2 + mb h - ib2 + mb 

+ k2- qb2Y - mf' {k2 - qb2Y - mf' ' 



xp(g) 



(gel + ge2) 



2 75 \ Vs'{qc2) 



in' \- ( \- f / h + h-ib2 + mb . h-ib2 + mb . 



xp(g) 



(gel + ge2) 



2 75 \ Vs'{qc2] 



fn \- ( \ - I / j\i ibi -h + ^b h-4b2 + mb 
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FIG. 2: Feynman diagrams that can be directly grouped into the bb subset. Here i and j are the 
color indices of c and b respectively. 



2/ 



Xpil) 



{QcI + qc2) 



2ls\'"s'iqc2) , 



M2d = {C3ij)us{qbiyi 



''bi - It2 + mb Iti- + rrib 



(2^ V^^ {qbi-hy-mf'ih-qb^y 



m 



2" 



xp{q) 



t75 Us'(?c2) , 



M2e = 



(gel + qc2Y ) 

Xp(?) 



{Cuj)us{qhi)i 



2ls\'"s'{.qc2) , 



d^q jjx^ jbi - ^■-> + Ai ~ i^'i ~ + "'-fe 



mt 



lis- 



xp{q) 



2l6\vs'{qc2) 



M2, = 



(gel + gc2) 



ib2 + "^6 la^l],^2l 



. (27r)4 \ (A;i + A;2 - 562)^ - (A;i + A;2)2 



{qci + gc2)' 



75 }Vs'{qc2) , 



Moh = 



{C2ij - Cuj)Us{qbi)i / TT^i , \ ^, n. , T"^^ 1^75 • 



(27r)4\(A;i + A;2)2 (^1 + ^2-561)2 



mt 



(g ^+ j2 7'5((^l - A;2)a5'Mi/ + (2/^2 + ki)ix9ua + (-2/?! - /C2)5'/.a) |fs'(gc2) ■ (5) 
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(0 



FIG. 3: Feynman diagrams that can be directly grouped into the cb or be subsets, where the first 
four diagrams belong to the cb subset, while the last four belong to the be subset. Here i and j 
are the color indices of c and b respectively. 

The third group, Fig|21 

ri - r \j\2 ibi-h + ^b .f d^q Xp{(1) h-ic2 + m^ 

lsVs'{qc2) , 

n - I ^- f d^q h-^b2 + mb Xp{q) h - ic2 + rric 



3d 



3/ 



^l'Vs'{qc2) , 

r' - { V f h - ib2 + mb Xpjq) jx, id -h+^c 

l5Vs'{qc2) , 

,Ai ibi -h + ^b . f d'^q Xpjq) h - ic2 + 



^ci - h + m^ 



^ wAi ibi-h + ^b .[ d^q xp{q) 



{q,,-k^y-ml J {2n)^ '\h - q^i - qc2r' {qa ~ k^f - ml 

l5Vs'{qc2) , 

n - I \- [ d'^q h-ib2 + mb . Xp{q) h - ic2 + 

M3, = C,,uM.J ft-fa-to)^^'ft-faP-mr 

^2^Vs'{qc2) , 

M r' (n V- [ J-!L^, h-ib2+mb Xpjq) jX2 ici-h+^c 

Ms. = C,,uMb.yj J^.^^^k,-qb,Y-ml^- {k. - q. - q.2Y^' {q.. - - ml 
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ntfrA — * — mtfSrl — * — rrrnit — * — rnwt — * — 



IP k 



mmnrB <- 

(a) 




qc2.i k, 




kiV n k 

TfTTfTit * 



mmnrrt — >- 

(e) 




qc2.i k, 



qc2.i ks 



(b) 



(c) 



(d) 




q-' ^^^^^ — 

t — ^-^ rrSri — m«3t^ — 



(f) 



-> mmrt — > 

qbi.j k, qt„,j k 

(g) 




qb,,j 



(h) 



FIG. 4: Feynman diagrams with only one three-gluon vertex, which can not be directly grouped 
into the cc, bb, cb and be subsets. Here i and j are the color indices of c and b respectively. 



The fourth group, FigEl 



(6) 



4a 



4b 



4c 



4d 



4e 



4/ 



(C4ij - C2ij)Us{qhl)i / 777^75 



{ki - 2gbi - 2qh2)^,ga5 + (Qbi + qb2 - '2ki)s9f,a 

74 



Xp{q) A2 ^cl 

[75: ^ 



1^2 + 



[{h + qbi + qb2)a9,i5 + 

Vs'iqc2) , 



( C75.,)^.(g,0^ y ^2vr)4 " (q,, + q..^' (q,, - k^f - ml (k, - q,, - q,,Y 
{{ki + qbi + qb2)a9tiS + {ki - 2qbi - 2qb2)f,gaS + (%i + qb2 - 2ki)sgf,a)vs'{qc2) , 

d^q xp{q) 



{Csij - Cuj)us{qbi)i 



[is: 



(27r)4 (g5i + qi,2y {k2 - qbi - qb2Y 

{k2 - 2qbi - 2qb2)f,ga5 + (qbi + %2 - 2/C2)5fi'^a 



k2 + qbi + qb2)a9iMS + 



i-qc2 + rric jx, I . 

l^Vs'{qc2) 



{ki - qc2Y - ml 



d^q xp{q) 



+ mr 



( C,.-,)uMbi)r J ^2^-^^75 (^^^ ^ ^^^^,fi ^^^^ _ _ _ _ 

{{k2 + qbi + qb2)a9t,s + {k2 - 2qbi - 2qb2)f,gaS + {qbi + qb2 - 2k2)sgf,a)vs'{qc2) 

xp{q) 



7aeiM 



fr^ n \- ( \- f ^bi -h + ^b 

(CH,-c3.,K(gbi>y j^j2 ^q^^ - k2f - mijk; 

{{ki + qci + qc2)agt,s + {ki - 2gci - 2gc2)/.fi'a5 + (gel + qc2 - 2ki)sg^a)vs'{q, 



qci - qc2y (gel + gc2)^^'^ 



{C5ij)Us{qbi)i 



d^q 



h-^b2 + mb . xp{q) 



r75 



(27r)4 {ki - qci - qc2Y {k2 - ^62)^ - ml ^ (g^i + gc2)^ 
ki + qci + qc2)agtiS + {ki - 2qci - 2qc2),,ga5 + {qci + qc2 - 2ki)sg^o)vs'{qc2) 
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FIG. 5: Feynman diagrams with two three-gluon vertices or with a four-gluon vertex, which can 
not be directly grouped into the cc, bb, cb and be subsets. Here i and j are the color indices of c 
and b respectively. 

Mig = {C2ij - C4ij)Us{qbi)i / JTT-Uh -f rT-2 271 ^7 T • 

{{h + Qci + qc2)agt,s + {h - 2gci - '2qc2)tiga5 + (gel + qc2 - '2k2)5g^a)vs'{qc2) , 

J (27r)4 {k2 - qci - qc2r ih - Qb2r - [Qci + qc2r 



{k2 + qcl + qc2)agtiS + {h - 2gcl - '2qc2)tiga5 + (gel + gc2 - 2fc2)<5fi'Ma)^s'(gc2) • (7) 

The fifth group FiglH 

M - rr^ n n \- i \ - [ J^l^^ Xp(g) lp^2l 

J (27r)4 {qf,^ + qf,2y {k2 - gcl - ge2)^ (gel + ge2)^ 

2A;i - qbi - qb2)ag,i5 + (2gw + 2gb2 - h)t,ga5 + (-gw - gfe2 - ki)sg,ia] 



[h + qcl + qc2)5gv(} + {k2 - 2gci - 2qc2),ygf35 + {-2k2 + gei + qc2)i3gsu)vs'{qc2) , 

in n n \- f \- f ^'^^ ^°^2^ ^^(^) ^^^^'^ 
M^b — \yuj — (^'iij—(^bij)'Us\qbi)'i I 777^7 \ T^Ti v)l \ V) ' 

J (2vr)4 (g;,! + q^^2f {h - gel - ge2)2 (gel + ge2)2 
{{2k2 - qbl - qb2)agfi5 + (2gw + 2gb2 - k2)^,gaS + (-gbl - g62 - k2)5gfia) ■ 

{{kl + qcl + qc2)sgvp + {ki - 2gcl - 2qc2)ugf55 + (-2/^1 + gel + ge2)/3fl'<5i.)^^s'(gc2) , 

M,. - (Cn, + C,., - C3,- - C2.,)uM.l)^j ^^q,,+q,,y^k, + k2r{qci + qc2f ' 

+ k2 + qbl + qb2)pgas + (^1 + k2- 2qbi - 2qb2)sgai3 + {qbi + qb2 - 2ki - 2k2)agi3s) 
{{kl - k2)5gut, + (2A;2 + ki)f,g^s + {-2ki - k2)ugs,})Vs'{qc2) , 
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((^^21^ + Csij Clij — Ciij){g^j3guoi — gfiagufs) + {C^ij — Ciij + Csij) ■ 

{9liu9l3a — 9tia9up) + (C'sij " C'2jj + C iij){g ^j^ygafj — giil39va)]Vs'{(lc2) ■ (8) 



In Eqs. (4-8), /ci, and ei, €2 are the momenta and the polarization vectors of the gluons; 
gel, Qbi are the momenta of c and h quarks and gc2, '?62 are the momenta of c and 6 anti- 
quarks, respectively. Note that in all the above equations we have omitted the factor g^ 
(the fourth power of the QCD coupling constant), so we should consider it when evaluating 
the final result. For convenience, we group the terms (Feynman diagrams) according to the 
character of the gluon attachment to the fermion lines. The details will be explained in the 
next sections. 

Under the non-relativistic approximation, for the weak binding system of (c6) we have 
MB^ = MBi=M^m, + m,, q^, = , = _p ^ (9) 
and the wave function Xp(Q') can be written as 

xp(g) = <^(^)^(«75 + m-sMP + M) , (10) 



where a = 1, /5 = for the pseudoscalar Bc{\I^Sq\) and a = 0, /3 = 1 for the vector 
-B*([l^S'i]). The radial part of the momentum space wave function 0(g) is related to the 
space-time wave function at origin by the integration: 

B. Motivation and basic formulae for decomposing the gluon self-coupling vertices 



In Ref. 



21[, a method is proposed for treating the amplitude of the processes, which 



contain massless fermions and non-abelian gauge boson(s), in order to make the result 
compact and to avoid numerical cancellations between very large numbers. The authors of 
Ref. ^ group the Feynman diagrams of the concerned process into gauge-invariant subsets 
according to how the lines of the gauge bosons attach to the fermion lines. They then choose 
convenient gauges for the subsets independently of each other (not a unique gauge for the 
whole process), and finally they obtain a compact result, when all terms of the amplitude are 
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written according to the helicities of the fermions. The key point of the approach (hehcity 
technique) is that when all the massless fermions are written on helicity state basis, then it is 
straightforward to choose convenient gauges for all subsets. If massive fermions are involved 
in the concerned process and one wishes to use the same technique, then one needs to 
generalize it. Some rearrangements, such as replacing one massive fermion by two massless 
fermions etc, should be made, and the gauge choice for each subset is more complicated 
than in the massless cases. In the present case for the subprocess gg ^ Bc + b + c, a. unique 
gauge for the whole amplitude is more practical. Apart from the gauge choice, we apply the 



techniques of Ref. j2l| as much as possible to make the program compact. The polarization 
vectors of the gluons are replaced by 7-matrix elements of massless fermion helicity states and 
the color factors are dealt with independently from the Dirac 7-matrix strings. Furthermore, 
we try to make the amplitude more symmetric by applying a decomposition of the terms 
when writing the program. This is achieved by decomposing first the terms, which contain 
three- or/and four-gluon vertices, into terms without self-interactions of gluons. Then, 
according to the structure of the contained fermion lines, some of the decomposed terms 
are chosen as 'coordinators', which would be called as the typical ones, so that all the 
other terms, referring to the Feynman diagrams, may be 'expressed' by the 'typical ones'. 
Therefore when writing the program, only different 'typical Feynman diagrams' need to 
be written precisely, while the non-typical ones may be generated by means of the typical 
ones according to the relationship (expression). In this section, we show how to decompose 
the terms in the Feynman diagrams containing three- or four-gluon vertex(ices) into terms 
similar to the typical ones, disregarding the difference in color factors, and list the results 
for the typical Feynman diagrams. 

First of all, let us introduce a massless fermion with an arbitrary light-like reference 
momentum q {q^ = 0) and its relevant helicity spinors {\q±) = ^^^l*?) and (^\q) = 0), and 
then construct the requested massive fermion four-spinors u{p) and v{p) with momentum 
p {p^ = m?) in terms of as follows: 

Us{p) = 2- {i) + m)\qh) , 

^ zp ■ q 

Vsip) = -^{^ - m)\q.,) , (11) 

where s = ±| is the spin of the massive spinors, while h = ± is the helicity of the massless 
spinor. For convenience we adopt the explicit form for the polarization vectors of the gluon 
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with momentum k as in CALKUL. When the hehcity states \k±) and of two massless 
fermions are defined as \q±), but the fermions have a hght-hke momentum k [k"^ = 0) and a 
referred one q' (g'^ = 0) respectively, the polarization vectors of the gluon with momentum 
k may be represented as follows: 

W-\l,\k-) 



^uik,q') 



V2{q' ■ k) ' 
V2{q' ■ k)* ' 



r{k.q) = j/^i\k+){l'^\ + W-){k-\). (12) 

Throughout the paper {q' ■ k) and (g' ■ k)* denotes and its complex conjugation 

respectively. 

Because the light-like momenta q and q' can be chosen arbitrarily, we choose them to 
be the same as a light- like momentum go 5 which is the reference momentum for all massive 
spinors and gluon polarization vectors. We take one gauge for the whole set of 36 Feynman 
diagrams, that is convenient here and different from CALKUL where different gauges are 
taken for different gauge-invariant subsets. 

The gluon self-coupling vertices do not adapt well for being simplified by using the po- 
larization vector in calculating the amplitude. Hence an effort to replace the part with 
the gluon self-coupling of the diagrams by the so-called QED-like ones (without gluon self- 
coupling) is described here. This approach is not straightforward, due to the fact that the 
gluon self-coupling diagrams with massive quarks can not be mapped completely into the 
QED-like diagrams as in the case of the massless quark condition|2l|. In the massive quark 
case, some additional functions need to be introduced. Fortunately, for the 'simple' subpro- 
cess, these 'extra' functions for the three-gluon vertices are just parts of the diagrams with 
four-gluon vertices. Thus the diagrams involving four-gluon coupling vertices just need to 
be decomposed according to their color factors and there is no need to reduce their 7 matrix 
and spinor factors any further. 

To decompose Feynman diagrams with self-interactions into those without any, we need 
to deal with some of the so-called 'basic structures'. The number of the 'basic structures' 
required for a specific process increases with the number of the gluons involved in the 
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concerned process. For the subprocess, gg Bc{B*) + 6 + c, to convert the necessary parts 
to the 'basic structures', only the parts containing a three-gluon couphng vertex have to be 
decomposed. Thus let us outline the decomposition for the concerned subprocess. 

The decomposition of a three-gluon coupling vertex is shown in Fig|Ul(the first structure): 
a three-gluon vertex /""'"^T^usi^i, K, ^2) through a gluon propagator —ig^^i/K"^ couples to a 
quark line with a quark-gluon-quark vertex T^-^^i, where K = —{ki + ^2) = —{Q + Q'), and 
jafec^ Tfe are color factors at the two vertices. It is 

M;iih, k,, Q, Q') = -g',r''nT,,sik„ K, k,)^^^^^,, , (13) 

where Q and Q' are the momenta of the fermion legs at the vertex. Since 

T^^s{ki, K, = {ki - K)sg^„ + {K - k2)^gu5 + {k2 - ki)^gf,5 , (14) 

so the precise expression is 



M^^K^i, ^2, Q, Q') = glr^'T,— [{k, + Q + Q')si, + k.^^s - - iMis 

((^'-m)(757^-^^5) + 

(7m7<5 - g^l&){Q + m) + k25'Jf, - fci^7<5^ • (15) 

If the factor g^, the propagator scalar factor and the color factor are disregarded, and the 
symbol ^ is used for an equality modulo these factors, we have 

+ (7m75 - fl'M<5)(^ + m) + A;257m " ^1m7<5 H , (16) 

where m is the mass of the fermion. Here the gluons and fermions may not be on-shell. If 
the first and the second terms on the right-hand side are embedded in the diagrams, the 
basic QED-like diagrams may be obtained, while the rest terms ■ ■ ■ on the right-hand side 
can be absorbed into several extra functions. 

To relate to the 'basic structures', two further decompositions are shown in FigI7| 
(the second structure) and Fig|Hl (the third structure) respectively: a three-gluon vertex 
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+ 



FIG. 6: The three-gluon coupling vertex is decomposed as in Eq. (|16|) : the first two terms are the 
'basic QED-hke' terms and the 'remaining' terms are expressed by several extra basic functions. 



FIG. 7: Reduction of the basic structure with a three gluon vertex: the first two terms correspond 
to the basic QED-like diagrams and the symbol 'x' means a quark gluon vertex T^^a. 

fabcTf_iu6{ki, ki + k2,k2) is coupled to a quark line, which has two vertices of quark-gluon- 
quark Tb'-fi, and T^'ja- To be precise, we label the latter vertex with the symbol 'x' on the 
quark line. In the present case the two external quark lines are both on-shell and then these 
two structures can be simplified by using the on-shell conditions 




+ c1 .X 



m) = 0, {g + m)v{Q) = 0. 



(17) 
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The contribution from the corresponding part in FigI7| is: 



gsr^'nTMQ'H 



a 



(fci + k2- QY - w? 



,vv' 



T^ys{ki,K, k2) 



X2 



lu'v{Q) . 



(18) 



With Eq.(|16p and the on-shell condition Eq. (jl7|) for the two external fermion legs, we 
obtain 



MX^ki, Q, Q') - uiQ'hcih + + {isih - ^ + "^)7m - + m)^s)v{Q) 

+ {ml + ml + 2ki ■ k2 - 2Q ■ ki - 2Q ■ k'^u{Q')'^a{l5liJi - gt^s)v{Q) + 



where '~' means that the factor g^, the propagator scalar factor and the color factor have 
been omitted. The first term is for the basic QED-like diagrams. The second term (cl ■ X) 
where 



will be treated below. The third term in fact, does not contribute to the amplitude, no 
matter whether the gluons are virtual or real. The proof is that when both the concerned 
gluons are real, it is easy to show that the remaining terms give zero contribution by using the 
relation e{ki) -ki = (Z = 1, 2). When one of the concerned gluons is virtual, the gluon with 
momentum k2 for example, k25 will always couple to a 'simple' quark line as u{R)'ysv{R') at 
the lowest order for the amplitude of the considered process gg Bc{B*) + c + b, where u 
and V are the quark and anti-quark spinors, and k2 = R + R' (the momenta R, R! satisfy the 
on-shell condition: R^ = R'^ = m? , where m is the mass of the quark). It is easy to show 
that the contribution of the remaining terms is zero by using a similar on-shell condition as 
Eq.()17|). Therefore, in Figdand in the following FiglSl the remaining terms have not been 
shown explicitly. 



+ (cl-X) + --- , 



(19) 



cl = m\ + ml + 2ki ■ k2 - 2Q ■ ki - 2Q ■ k2 



and 
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c,6 



c,(5,k2 




+ c2.Y 



FIG. 8: Dividing the basic topology including the three gluon vertex, where the first two terms 
correspond to the 'basic QED-like' diagrams and the symbol 'x' means a quark gluon vertex T^^a. 

In the same way as above, the contribution from the corresponding part in FiglHlis: 



MXih, k2, Q, Q') ^ u{Q') {isi^' -h + Mh, - -h + Mhs) -h~h + M)- 

lav{Q) + {rn\ + ml + 2ki ■ k2 - 2Q' ■ ki - 2Q' ■ k2)u{Q'){-f5lt^ - Qt^s) ■ 
'javiQ) + u{Q')(-k25i^{h + h) + Klsih + h) + 2A;25Q; + 



where again the factor g1, the scalar factor of the propagators and the color factor have 
been omitted. The first two terms correspond to the basic QED-like diagrams, and the term 
which is expressed by (c2 ■ Y) is defined as 



Similarly as for the second structure in Eq.(19), the remaining terms in the present structure 
contribute nothing, thus in FiglHl they have not been shown explicitly. 

Having made all the preparations above, and performing all the possible interchanges, 
such as gluon exchange, quark and anti-quark exchange, we relate each of the terms in 
the diagrams to a combination of those with the basic QED-like diagrams ('coordinators'), 
including the introduced two 'extra' functions as well. 




(20) 



c2 = ml + ml + 2ki ■ fca - 2Q' ■ ki - 2Q' ■ k2 , 



where F is a new extra function 



Y = m(Q'){7<57/. - gf^shaviQ)- 
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C. The decomposition 



The Be meson is a double-heavy weak-binding state. According to pQCD each term of 
the amphtude for the subprocess may be factorized into two factors: that of perturbative 
gg ^ b + b + c + c (all the quarks are on shell) and that of non-perturbative c + 6 — Be- 
The binding wave function in the Bethe-Salpeter framework may be used to dictate the 
non-perturbative one, and approximately written as Eq. lfTm) . To carry out the factorization 
of the amplitude, one may apply Eq.lfTIH) and the two equations: 

^ci + mc = J2 "(^ci, s)u{qci, s) ~ M) , 

s 

ib2 = v{qb2, s)v{qb2, s) ~ -a2(f + M) (21) 

s 

to each term in the amplitude in Eqs.(@]|Hl) with the corresponding factors Xp(Q')- 
Then the general structure of the amplitude in 'explicit helicity' form turns to 

^(A„A4,A5,A6)^^^^^ gb2, gel, qc2, K h) = E aX,Z}i4"/''^'"^'"^'"^'"'')(gw, g,2, qa, qc2, h, h) ■ 

A2,A3 

^24':(sf)(gh2,gci), (22) 

where i = 1, ■ ■ ■ , 36 (or labelled as Feynman diagrams: la, 16, ■ ■ ■ , 5d, respectively), (j = 
1, ■ ■ ■ , 6) denote the helicities (spins) of the quarks and gluons respectively appearing in the 
two factorized 'processes' gg b + b + c + c and c + b ^ Be- Note that from now on, 
we change the notation of the helicities of the particles in the processes as: Ai denotes the 
helicity of b, A2 that of b, A3 that of c, A4 that of c; whereas A5, instead of Ai in Eqs.(4-8), 
denotes that of gluon-1 and Ae, instead of A2 in Eqs.(4-8), denotes that of gluon-2. Here 
Ci, Xi denote the color factor and the scalar factor from all the propagators as a whole for 
the zth-diagram, respectively. B''pl'^^'^'^'^^'^'"^'^\qbi,qb2,qci,qc2,ki,k2) and 5^^^''^j)(gb2, gci) 
are the amplitudes corresponding to the 'free quark part' g{ki, X^)g{k2, X%) b{qbi,Xi) + 
b{qb2, ^2) + c(gci5 A3) + c(gc25 A4) (all the quarks are on-shell) and the 'bound state part' 
c{qci, A3) + b{qb2, A2) Bc{B*), respectively. Substituting Eq. (fTT|) . we have 

B^pl'^^'^'''^'''^'"^''\qbl,qb2,qcl,qc2,h,k2) = (goAi |(^W + "^fe)rii(^b2 - "^fe)|«?0A2) • 

(goAaK^ci + mc)T2i{ic2 - mc)\qQ\^) , (23) 

<(B2!(?^2,gci) = (goA2l(«75 + M^.))^v7^^(0)l^o^3), 

(a = l, /3 = 0, for B,[l^S,]- a = 0, /3 = 1 , for ^^[l^Si]) (24) 
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where \qoK) and {qoKl {K = ±, r = l,---,4) are the introduced hehcity states of the 
massless fermion go which specifically relate to those of massive fermions with the mo- 
mentum p and mass m as in Eq.(|lip. rij2j are the explicit strings of Dirac 7 matrices 
corresponding to i-th Feynman diagram which contain the gluon helicities A5 and Ag. 

Di = -y= 7= 7= 7= — and Do = 7^ — are the two common normal- 

ization factors. 

The function D2B^^^^^l}^{qi,2, Qci) , which contains the bound state wave function, is: 



D2Bl^'''\q,2, gel) = ^^^1^ ^A2.A3(^A2- - 5a2+) (25) 

for and 

[Qb2,qcl) - o / ^A2,A3l^A2+ + 0A2-j 5 + 



2^mbmc ' \ P-qo J 2^mbmc \2P ■ qo) 

(goA2l/(s.)f|goA3) (26) 

for Bl[^Si]. 

These 36 functions, B'^pl'^'^'^^'^'^'^^'^^\qbi,qb2-,qci-,qc2iki,k2)., can be constructed by nine 
'basic ones' which correspond to FigOt', FigOb, FigOfc, Fig Hi, FigOl (two basic functions), 
Fig^, FigO^ and Fig^, and by performing all the possible interchanges of the initial 
gluons and the two quark lines. In fact, the functions which correspond to FiglHfc and Fig Oil 
can be obtained from those in FigOti' and FigOb by interchanging the initial gluons and the 
quark lines in the diagrams, and in the following section we will show that the functions 
corresponding to Fig|T^ and FigHJ:; can also be expressed by other seven 'basic functions', 
although here we still treat those four functions as 'basic ones', in order to treat them on 
an equal footing. 

We use Em,j,k{qbi, qb2, qd, qc2, ^1,^:2), (m = 1, 2, ■ ■ ■ , 9; j = 1, ■ ■ ■ , 4; k = 1,2,--- 64) to 
denote the 'basic functions', where k denotes 64 possible helicities (spins) corresponding 
to possible 'values' of (Ai, A2, A3, A4, A5, Ag) as shown in Table HI and the correspondences 
of the functions to the Feynman diagrams are shown in Table |n] and Table IIIII Here j 
means the four possible interchanges: 1 means identical (without any interchange), 2 means 
interchange of the gluons, 3 means interchange of the quark and the anti-quark, and 4 means 
interchange of the gluons and the quark and anti-quark. When applying the interchange 
among the particles, the symmetries of the decomposed amplitude guarantee that every 
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TABLE I: The correspondence between k = 1, • • • , 64 and Ai = it, A2 = ±, A3 = it, A4 = it, A5 = 



lb, Ag = ±, which stand for the hehcities of the particles in the process. 



k 


Ai 


A2 


A3 


A4 


A5 




k 


Ai 


A2 


A3 


A4 


A5 


As 


k 


Ai 


A2 


As 


A4 


A5 


As 


k 


Ai 


A2 


A3 


A4 


As 


Ae 


1 


+ 


+ 


+ 


+ 


+ 


+ 


17 














33 


+ 


+ 




+ 


+ 


+ 


49 






+ 








2 


+ 


+ 


+ 


+ 


+ 




18 












+ 


34 


+ 


+ 




+ 


+ 




50 






+ 






+ 


3 


+ 


+ 


+ 


+ 




+ 


19 










+ 




35 


+ 


+ 




+ 




+ 


51 






+ 




+ 




4 


+ 


+ 


+ 


+ 






20 










+ 


+ 


36 


+ 


+ 




+ 






52 






+ 




+ 


+ 


5 


+ 


+ 


+ 




+ 


+ 


21 








+ 






37 


+ 


+ 






+ 


+ 


53 






+ 


+ 






6 


+ 


+ 


+ 




+ 




22 








+ 




+ 


38 


+ 


+ 






+ 




54 






+ 


+ 




+ 


7 


+ 


+ 


+ 






+ 


23 








+ 


+ 




39 


+ 


+ 








+ 


55 






+ 


+ 


+ 




8 


+ 


+ 


+ 








24 








+ 


+ 


+ 


40 


+ 


+ 










56 






+ 


+ 


+ 


+ 


9 


+ 






+ 


+ 


+ 


25 




+ 


+ 








41 


+ 




+ 


+ 


+ 


+ 


57 




+ 










10 


+ 






+ 


+ 




26 




+ 


+ 






+ 


42 


+ 




+ 


+ 


+ 




58 




+ 








+ 


11 


+ 






+ 




+ 


27 




+ 


+ 




+ 




43 


+ 




+ 


+ 




+ 


59 




+ 






+ 




12 


+ 






+ 






28 




+ 


+ 




+ 


+ 


44 


+ 




+ 


+ 






60 




+ 






+ 


+ 


13 


+ 








+ 


+ 


29 




+ 


+ 


+ 






45 


+ 




+ 




+ 


+ 


61 




+ 




+ 






14 


+ 








+ 




30 




+ 


+ 


+ 




+ 


46 


+ 




+ 




+ 




62 




+ 




+ 




+ 


15 


+ 










+ 


31 




+ 


+ 


+ 


+ 




47 


+ 




+ 






+ 


63 




+ 




+ 


+ 




16 


+ 












32 




+ 


+ 


+ 


+ 


+ 


48 


+ 




+ 








64 




+ 




+ 


+ 


+ 



term may be expressed by the basic functions with helicities and momenta of the particles 
in the process. For all the functions Em,j,k, with m,j being fixed, they are related to each 
other by proper complex conjugation with or without changing the whole sign. To write the 
program, and to apply these relations of the element functions among different helicity states 
conveniently, we take the correspondence between k and (Ai, A2, A3, A4, A5, Ae) as described 
in Table m Thus we have now 

^{Ai,A2,A3,A4,A5,A6)j^^^^^ gfe2, Qcl, qc2, h, h) = B^Fi iQbl, ^62, Qcl, qc2, h, 

9 4 

~ X/ X/ fi,m,jEm,j,k ■ (27) 
m=l j=l 

Here i and k denote the 36 Feynman diagrams (terms) and the 64 possible helicities of all the 
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TABLE II: The expansion coefficients fi,m,j for the functions Bp- {qf)i,qb2, qci,Qc2, ^i, ^2) which are 
grouped into the cb subset directly or indirectly through a proper decomposition (the coefficients 
fi,m,j are not listed here if they are equal to zero in a whole row). 
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-2(9c2 • ^1) 


fid,m,j 





-1 
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2(gci • ki) 




1 


-1 








2(^61 • k2) 


fig,m,j 








1 


-1 


-2{qb2 ■ k2) 


f5b,m,j 


1 


-1 


-1 


1 


Sb + Sc- S2 



particles in the two factorized processes. The correspondence of hehcities on k = 1, ■ ■ ■ , 64 is 
given as described in Table lU The 36 functions i = 1, ■ ■ ■ , 36 are labelled as the 36 Feynman 
diagrams, i.e. i = la, 16, ■ ■ ■ , 5d. The coefficient fi^mj corresponding to QED-like diagrams 
can be directly read out easily, but as for the diagrams involving three- and four- gluon 
vertices we decompose them by applying the formulae obtained in the subsection B. 

D. The amplitude B'-p]{qbi,qb2,qci,qc2, ^2) 

All the terms for the subprocess may be divided into four subsets according to the manner 
how the gluons attach to the fermion lines. For the cc set, both gluons attach directly to 
the c-quark line; for the cb set, the gluon-1 attaches to the c-quark line while the gluon-2 
attaches to the 6-quark line; for the be set, the gluon-1 attaches to the 6-quark line while 
the gluon-2 attaches to the c-quark line; for the bb set, both gluons attach directly to the 
6-quark line in the relevant Feynman diagrams. The diagrams involving three- and four- 
gluon vertex (vertices) which are shown in FigEl and FiglHl are decomposed by applying 
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the decomposition formulae in subsection B. Hence each of the terms corresponding to 
the diagrams in general turns to four terms and then they may be put into four subsets 
separately according to the resulting characteristics of the decomposed term. The results 
are shown in Table O and Table 11111 In Table O and Table 11111 we list the exact results for 
the cb and cc sets, and the results for the be and bb sets by performing interchanges of the 
initial gluons and the final quark lines. In Table ITTl and Table UTTl one may see, for example, 
how the Feynman diagrams in FigsEJ:;, EJi and FigEb are decomposed, and how each of the 
decomposed term is divided into a c6 or a cc set. We summarize the relations among the 
functions 5^ (g,i 

! 962, Qci, qc2, ki, ^2) and show them in Table HIl and Table Hill respectively. 
In Table HTl the functions Bp] {qi,i, qh2, Qcu (lc2, ki, are hsted either directly or through a 
proper decomposition in order to relate them to the cb subset through functions Emj,k and 
suitable coefficients fi,m,j- In Table HTTl the functions Bp]{qbi,qb2,qci,(lc2,ki,k2) are listed 
either directly or indirectly through a proper decomposition relating to the cc subset via 
suitable coefficients fi,m,j- In addition, in Table ITTT| we have decomposed the matrix element 
of FigOi, as shown in Eq.(jHl), into three parts according to the three different color 
factors: C^di = C2ij + C'nj — Cuj — Ciij, C^d2 = C^ij ~ Cuj + C^ij, C^ds = C^ij — C2ij + C/^ij. 
The parameters used in Table HT] and Table IIIII are as follows: si, = (gti + qb2Y, Sc = 
(gel + gc2)^ S2 = {k2 - qbi - g62)^ A = 2(gci - qc2) ■ k2 + Sb, f2 = 2(gci - qc2) ■ ki + Sb, 
/s = 2ki ■ k2 - 2qc2 ■ ki - 2q^2 ■ ^2, h = 2A;i ■ k2 - 2q^i ■ ki - 2q^i ■ /cg- 

E. The 'basic functions' Em,j,kiqbi,qb2,qci,qc2,ki,k2) 

To obtain the nine 'basic functions' Em,j,k = 1, ■ ■ ■ , 9), let us first define the functions 
which correspond to various kinds of quark lines (different 7 structures of the fermion lines), 
where q\ = q2 = ""^^ and q^ = k'^ = ki = k^ = 0: 



/o(gi,g2, Ai, A2) 


= {qoxi 


/i(gi,g2,fc, Ai, A2, A3) 


= {qoxi 


/2(gi,g2,fc, Ai, A2, A3) 


= (goAi 


/3(gi,g2, A;, Ai, A2, A3) 


= {qoxi 


fi{qii q2, ki, k2, Ai, A2, A3, A4) 


= {qoxi 




ih- 



1(^1 +m)75(^2 -m)|goA2) , 

1(^1 + m)75(^ -<i2 + m)f^{k, go) ((^2 - m)\qox2) , 
1(^1 + m')^^''{k,qo){(ii + m)-fsi<i2 - m)\qox2) , 
1(^1 + m)/^»(A;,go)(^2 -"^)|goA2) , 
1(^1 + m)-fs{h + h-i2 + m)f-'{ki, go) • 
- i2 + m)f'^{k2,qo){i2 -m)\qQx^) , 
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TABLE III: The expansion coefficients fi,m,j for the functions Bp- {qhi,qi,2, qci,qc2, ^1,^2) which are 
grouped into the cc subset directly or indirectly through a proper decomposition (the coefficients 
fi,m,j, which are equal to zero in a whole row, are not listed here). 
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f5{qi,q2,ki,k2,Xi,X2,X3,X4) = (goAiK^i '(A:i,go)(^i - h +m)js ■ 

{h -42 + m)/^*(A;2, go) (^2 - "^jkoAa) , 
/6(?l,?2,^l, A;2, Ai, A2, A3, A4) = (goAiK^i +H^^^(^i>?o)(^i - 1^1 +"^) • 

i^*{k2:qo){ii - h- h + i^)i&{h - i^)\qo\2) ) 
/7(gi, ^2, A;2, Ai, A2, A3, A4) = (goAiK^i + rn)^5f''{ki, qo)f*{k2: go)(^2 - m)\qox2) , 
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/8('?i,g2,fci,fc2,Ai,A2,A3,A4) = (goAiK^i +"^)/^'(fcl,go)/^*(fc2,go)7<5(^2 -m)\qox^) , 

/9(gi,g2,/i;i,/c2,Ai,A2,A3,A4) = (goAiK^i + 50)75/^^(^2, go)(^2 -"^)|goA2) • (28) 



There are several ways to deal with these fermion lines. In Ref. a systematic way for 
doing the massive case is proposed. Here we take another and more 'direct' approach. When 
the massive fermions have time-like momenta {i = 1, 2) and are directly connected to 
|goA,) or (goAil as in Eq. (j28|l . we may introduce the light-like momenta by defining 

(29) 



-qo- 



2qi ■ go 

Then for massive fermions can be replaced by the massless ones, ^/j, without any conse- 
quences. This is due to the relations 



^olgoA,) = 



(30) 



or 



(goAj^o = 0. (31) 

If the massive fermions with momenta {i = 1,2) are not directly connected to |goAj) or 
(goAj, but to another light-hke spinor, say |go;^.) or (goAj, then we may do the same thing 
if the momentum go and the spinor |goA,) or (goAj are replaced by q'^ and |goA,) or {qQXi\j 
accordingly. In this way we can turn the massive terms into mass, 
can be dealt with similarly as in the massless cases [ifll. ^ 
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ess terms, and then they 
2J. 



The nine basic functions may be expressed in terms of the above ten functions as 
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(32) 
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Applying the exchange symmetries between the two gluons and among the quarks, the 
following relations may be obtained: 

El^S^k = -^4,2,^, -E'2,3,fe = -E'2,2,fe, -E'3,3,fc = -£'3,2,*; , 

E4,3,k = -E'l,2,fc; -^5,3,*; = -^5,2, fe; -^1,4,^ = -^'4,1, A; , 

E2,4:,k = -S2,l,fe, -£'3,4,fc = -E'3,l,fc; -^4,4,^ = -£'l,l,fc , 

-E'5,4,fc = E^,l,k-> -^9,4,fc = -E9,l,fc + -£'9,2,fc " £-9,3,*: • (33) 

Furthermore, for the diagrams involving three gluon vertices, by using a proper decomposi- 
tion and the results in Tables HTl and IllH EQ j k may be replaced as 

-E'6,l,fe = -^7,1, fe + 2gc2 ■ ^2-£'9,l,fc ~ -^3,2, fc + -£-1,2, fe , 

-£'6,2, fc = -£'7,2,fc + 2gc2 ■ kiE() 2,k — -£'3,l,fe + -£'l,l,fe , 

-E'6,3,fc = -£'7,3,fc + 2g;,2 " ^2-£'9,3,fc — -£'3,4,^ + -£'l,4,fc , 

-£'6,4,fc = -E7,4,fe + 2gc2 ■ ^1-^9,4,^; — -^3,3,^ + -Ei^s^fc ; (34) 

and Ej^j^k may be replaced as 

-£'7,1,^ = —Ei^^k + -£'2,1,A; + -E8,l,fc — 2gcl " ^1 (2-^5,1, fe — 2E^^2,k + -^9,1,^) 5 

-£'7,2,fc = —Ei,2,k + -E'2,2,A; + -£'8,2,fc ~ '^Qcl ' ^2(2-^5^2, fc " 2i?5^i_fc + EQ 2,k) , 

E7;3,k = —Ei^'i^k + -£^2,3,^ + -£'8,3,fe — 2gfei ■ ki{2E^-i^k — 2E^^i^k + Eg^s,k) , 

-£'7,4,fc = — -^4,4,^; + -^2,4,^ + -E8,4,fc — 2gj,i ■ A;2(2i?5^4 — 2i?5^3 + £^9^4 fc) . (35) 

Thus the seven kinds of basic functions Emj,k ijn = 1, ■ ■ ■ , 5, 8, 9, j = l,---,4, k = 
1, • ■ ■ 64) may be written in a very compact form. As an explicit example, -Em,!,! {"^ = 
1, ■ ■ ■ , 5, 8, 9) is shown in one of the appendices. 

F. Rearrangement of the color factor for the amplitude 

As shown in section II, there are only five independent color factors, and we may choose 
them as Cmij {m = 1, ■ ■ ■ , 5), where i,j (1, 2, 3) are the indices of the final b and c quarks' 
colors respectively. Thus the whole amplitude may be rewritten as 

36 

j\^(Ai,A4,As,A6)(^^^^ 962, qa, qc2, fci, k2) = Yl M^^''^"^"^'\qki, qb2, qa, qc2, k,, fca) 

1=1 

5 

= Yl CmijM:i^"^''^"^'\qbi, qb2, qa, qc2, h, k2) , (36) 

m=l 
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where Cmij ("^ = 1, ■ ' ' > 5) are the five independent color factors, defined in subsection A. 
With Eq.(I221), M'^^^'^^'^'"^''\qbi,qb2,qci,qc2,ki,k2) (in short notation M^) can be further 
factorized as 



K = E C^^M^'''''''''''''''\qtU qb2, gel, gc2, h, fc2)I^2<(Bf)(^^2, gel) , (37) 



where M^^^'^^''^^''^'''^^''*'^''(gf,i, 5^2, gci, gc2, ^i, ^2) (in short notation M'p^) is the amphtude of 
the 2—^4 free quark process gg ^ c + c + b + b. 

Owing to the fact that each of the terms in the amphtude (^Fm) related to one 
of the 36 Feynman diagrams, the explicit formulae for M^{Mp^) {m = 1,- ■■,5) may be 
written down directly. With the nine kinds of basic functions Ei^rn,k Eqs. (l32llS^ . M'p^ may 
be written as 

A^i?! = — (2(X3a + + X55)£'i^l_fc — 2X4ei^'2,4,A: — 2X4c-E3,l,fc + 2X4e-E4,4,fc 





(38) 



M' 



F2 



2(X3e + X4a + Xr,a)Ei^2,k — 2X4gi?2,3,fc "~ 2X4ail^3^2,fc + 



2(Xia + Xig + X5c)-E6,l,fc ~ (Xlg + X5c)-E6,2,fc + X^c{Es^2,k — Eg^l^k) + 
(X2e + X2h)Es,3± — X2hE^A,k) — X5(i(2i?5_2,fe + -Eg.l.A; + -^9,2, fc) + 



2X4g(i?4 3 + 2E^ 3 kqbl ■ ki) + 4X4a£'5_2,fcgc2 ' ^2 + (Xig + Xsc) (-Eg^i^fc 
Eg,2,k)f3 + X5c(4i?5^2,fc — Eg^i^k + -£'9,2,fc)/4 + 4-E'5,l,fc ' {X^d — X^cfi) + 



X2h{4:E5^S,k — 4-E5,4,fc + -E9,3,fc — -E'9,4,A:)/6 — 2X5a{_E'2,2,fc + -£'3,2,fc 




(39) 



M' 



F3 



~^(~2(X4c + X5b)Ei^i^k + 2X4e-E2,4,fc + 2(X3c + X4c)-E3,l,fc + 
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-£'9,2,fe) — 2X4e(-E'4,4,fe + '^E^,i,k<lh\ ' ^2) — ^X/^cE^^i^k(lc2 ' ^1 

+X5c{Eg,l,k — Eg^2,k)f3 — i^lh + -'^5c)(4-E5,l,fc — 4-^5,2, fc + 

-E'9,l,fe — Eg^2,k)h ~ ^2g{Eg^3^k — Eg^^^kjfb + '^{{^3b + ^5b)E2,l,k + 

^56[^3,l,ik - E4,l,k - ^5,l,fe(-Sl + Sb + Se)]}) , (40) 

M^, = ^(-2(X4. + X5.)Ei,2,. + 2(X3, + X4.)£;3,2,. + 2X3.^4,2,.- 

2X5c£^6,l,A; + 2(X5c£^6,2,fe + (-^2o + -^2g)-£'6,3,fc " -'^2ff-E'6,4,fe + 
XicE^^l^k + X2cE7,3,k + {Xu + + X5c)-E'8,l,jt — (-'^l/i + -'^5c)-E'8,2,ifc) + 

-'^Sd (2-^5,2, fe + -Eg^i^fc + £'9,2, fc) + 2X4g(i?2,3,fe " -E'4,3,fc " 
2-E5^3_jtgbl • /Ci) — 4X4a£^5_2,A;?c2 " ^2 — X5c{E9,l,kE9,2,k) fs — 
— {^Ih + -'^Sc) (4-E'5,2,fe — -^9,1,?; + Eg^2,k)f4 + 4-E'5,l,fe(— ^s^f + 
(-^l/i + X^c)/^) + X2g{Eg^z,k — Eg^4^k)f5 + 2{(X3/ + X^a)E2,2,k + 
^5a[-E'3,2,ik - -E4,2,it - -E5,2,ik(-S2 + Sft + Sc)]}) , (41) 
^'f5 — -Dl(~(-'^56-E'l,l,fe) + (-^^4^ + -'^56)-E'2,l,fe " (-^^3^ + -'^4d)-E'4,l,fe " -'^3h-E'4,2,ifc " 
-'^2a(-£'6,3,fc + -E'6,4,fc) " ^le-E'8,l,fc " -^1/-E'8,2,A: " -^5d(-E'5,l,fc + -^5,2,fc 
—Eg^l^k — Eg^2,k) + -^4/i( — -E'1,3,A; + -E'3,3,fe " '^Er,^^ f,q^2 ' ^l) + 
X4/( — £'i^4^fe + -^3,4,*; — 2£'5^4^jfcgb2 • ^2) — "^XidE^^l^kQcl ' ^1 + 
X4b{E2Xk — E4^2,k — '^E^^Mcl ' ^2) " -^56(— -^3,1,/; + -^4,1,^ + 
-E'5,l,fe(— Si + Sb + Sc)) — Xr^a{El,2,k — E2,2,k — -E'3,2,it + 

^4,2,fe + £^5,2,fe ( -S2 + + Sc) }) , (42) 

where Xi are the scalar factors of the propagators, labelled by the corresponding Feynman 
diagram, si = {ki-qbi-qb2y, h = 2krk2-2qb2-ki-2qb2-k2 and /e = 2krk2-2qbrki-2qbrk2. 

G. Programme check for the subprocess gg Bc.{B*)+b+c and hadronic production 
of Be 

With all the above preparations, it is straightforward to write the programme for comput- 
ing the cross section of the subprocess. Applying the pQCD factorization theorem, various 
cross sections of the Sc(-B*) -hadronic production may be obtained by integrating the hadron 
structure functions over the cross section of the subprocess gg — > Bc{Bl) -\-h-\-c. 
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Various cross-sections of the subprocess gg Bc{B*) + b + c may be calculated, once 
the amplitude (with initial state color factors averaged and final state color factors summed 
up) has been computed. One only has to integrate over the proper phase space according 
to the requirements. Since the final state of the process is a massive three-body state, 
the phase space integrations are carried out numerically. As in the more general case, 
a Monte-Carlo simulation integration over the phase space is a practical solution, when 
one is averaging over the helicities of the initial gluons and spins of the quarks (and spins 
of -B* a gg B* + b + c is considered). To do the phase space integration, we first 
use the routine RAMBOS [26] to generate the requested phase space points (the energy- 
momentum conservation is guaranteed and some additional constraints for specific requests 
are matched). We then use the VEGAS H program (with necessary revisions to suit the 
present problem) to perform the integrations. The VEGAS program is useful for obtaining 
accurate total cross-sections, smooth distributions for the Pt|24] and the rapidity Y of the 
-Bc-nieson, and so on. When running VEGAS, the most important samples for the matrix 
element squared are taken first. Then by taking an adequate number of points for the 
integration, we obtain a final result, which is stable with respect to increasing the number 
of points and compatible with the requested statistical error. 

To check the program, we calculate the total cross-section of the subprocess, gg 
Bc{B^*^) + b + c, as well as the B^. and Y distributions for various subprocess center-of- 
mass energies, using the same parameters {e.g. a., rrib, rric, riiB^ and fBc{B*)) as in Refs.0,|8|. 
The results are compared with those in Refs.|2llSi and shown in Table HVl and Table El The 
Bc-pT and Bc-Y distributions for various subprocess center-of-mass energies are compared 
with Ref. j?! and shown in FigEl The integrated cross sections versus the center-of-mass 
energy of the subprocess for a fixed value of a.^ = 0.2 are also shown in FiglTUl One may see 
the agreement between our results and those in Ref . [A M from the tables and figures. Since 

pn 

the present programme and that used in Ref. [7|, |8|| are totally different, the agreement is a 
very solid check for both of them. 

According to pQCD, the production cross section is formulated 

dcr = J2j J dx2Flj^{xi, Hf) x Fj^^{x2, fJ.F)daij^B,x{xi, X2, fJ^p), (43) 

where F'^{x,^f) is this distribution function of the parton i in the hadron H, 
dcTij^Bcxixi, X2, fJ^p) is the cross section for the relevant inclusive production {i+j B^+X). 
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TABLE IV: Comparison of total cross sections for gg Bc{B*) + h + c with the corresponding 
results of Ref-i^]. The input parameters are rrif, = 4.9 GeV, nic = 1-5 GeV, Mb*c = rnb + "ic, 
fsc = 0.480 GeV, ag = 0.2. The number in parenthesis shows the Monte Carlo uncertainty in the 
last digit. The cross sections are expressed in nb. 





20 GeV 


30 GeV 


60 GeV 




0.6579(5) X 10-2 


0.9465(8) X 10-2 


0.7872(8) X 10-2 




0.661(7) X 10-2 


0.949(8) X 10-2 


0.782(9) X 10-2 




0.1606(1) X 10^1 


0.2460(3) X 10-1 


0.2033(2) X 10-1 




0.160(2) X W-^ 


0.244(3) X 10-1 


0.203(3) X 10-1 



TABLE V: Comparison of total cross sections for gg ^ B^. + b + c with the corresponding results 
of RefPl- The input parameters are = Sruc, Mb^ = 6.30Gey, fsc = 0.480GeF, = 0.2. The 
number in parenthesis shows the Monte Carlo uncertainty in the last digit. The cross sections are 
expressed in nb. 





20 GeV 


30 GeV 


60 GeV 


80 GeV 




0.6853(5) X 10-2 


0.9731(8) X 10-2 


0.7997(9) X 10-2 


0.6244(9) X 10-2 




0.686(2) X 10-2 


0.971(4) X 10-2 


0.793(5) X 10-2 


0.623(5) X 10-2 



With the pQCD factorization formula Eq. ()43|) . various cross sections for the hadronic pro- 
duction of the meson B(.{B*) can be computed by integrating over two more dimensions of 
the structure functions of the incoming hadrons. 

The programme for the hadronic production of the Bc{B*) mesons has also been checked 
by evaluating the hadronic production of B^. at Tevatron. The explicit example is to produce 
the Be meson by using the next to leading order running a^, the characteristic energy scale 
Q2 = s/4 (NQ2=1) of the production and the CTEQ3M set of parton distribution functions. 
The results for the distributions and the cross sections agree with those in Ref.j^. 
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III. THE PROGRAMME: BCVEGPY GENERATOR 



The programme BCVEGPY is a generator for hadronic production of Be mesons in form 
of a Fortran package, based on the dominant subprocess gg Bc{B*)+b+c. Concerning its 
implementation in PYTHIA, BCVEGPY is written in the same format as in PYTHIA (in- 
cluding common block variables). Thus it is easy to implement straightforward in PYTHIA 
as an external process, and in this way all the functions of PYTHIA can be utilized in 
connection with the use of BCVEGPY. 

A. Structure of the program 

The BCVEGPY Fortran package contains five files: bcvegpy.for, genevnt.for, sqamp.for, 
foursets.for and pythia6208.for. The routine bcvegpy.for is the main program of BCVEGPY 
which takes care of the necessary input parameters for the event generation and outputs a 
variety of results, such as the distribution of the transverse momentum and rapidity Y of 
the produced Be, etc.. In the routine genevnt.for, there is a function TOTFUN and five sub- 
routines: EVNTINIT, UPINIT, UPEVNT, BCPYTHIA and PHPOINT. The subroutines 
EVNTINIT and UPINIT are used for initializing the parameters when running BCVEGPY 
in the PYTHIA environment. The event generation starts by calling PYEVNT (a PYTHIA 
routine) just after calling the subroutines EVNTINIT and UPINIT. The routine UPEVNT 
(a PYTHIA user routine) is then called, that allows the implementation of an external pro- 
cesses. The functioning of the subroutine UPEVNT is determined by an input flag, which 
is a number read in from the file input.dat. When the input flag switches on, the generation 
of complete events may be carried out. UPEVNT calls the subroutine BCPYTHIA and the 
full information, i.e. the status code, the mother code, the color flow etc. for the final state 
particles {Bc{B*) meson and two jets b and c). The allowed phase-space for the production 
(i.e. energy-momentum consevation) controls the production completely. The phase-space is 
generated by calling the subroutine PHPOINT. The function TOTFUN computes the value 
of the integrand for the phase space integration by caUing the subroutine AMP2UP which 
is in sqamp.for. The subroutine AMP2UP is written according to the techniques described 
in the previous sections and its purpose is to compute the square of the amplitude for the 
hard subprocess gg — > Bc{B*) + b + c. In files sqamp.for and foursets.for, there are quite a 
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lot subroutines and functions, that are needed for calculating the square of the amplitude. 
All of them will be explained in Appendix C. When running the package BCVEGPY, the 
PYTHIA library must be linked; in particular, the generator is designed to be interfaced to 
the PYTHIA version 6.2 2^. For reference, we include the file pythia6208.for in our pack- 
age. In order to increase the phase space integration accuracy, one may apply the subroutine 
VEGAS (in sqamp.for) to optimize the sampling of the phase space points for phase space 
integration before calling PYTHIA. 



B. Use of the generator 

BCVEGPY can be used for generating a huge sample of Bc{B*) meson for hadronic 
collisions efficiently (based on the mechanism with gg —>■ Bc{B*) + 6 + c as the main sub- 
process). Furthermore, if the generator is implemented in PYTHIA, a complete event with 
two hadronized quark jets and a decayed Bc{B*) meson can be simulated. If one would 
like to simulate only the hadronic production of the Bc{B*) meson, BCVEGPY alone is 
sufficient. Furthermore, users may choose either the hadronic production or the subprocess 
gg Bc{B*) + b + c only, just by setting the value of the flag ISUBONLY equal to to 
switch on the integration for the parton distribution (structure) functions, or setting the 
flag equal to 1 for the subprocess only. Since the cross section of B^ production in hadron 
collisions is small compared to the B^ , B^ , Bg production, the efficiency for producing Be 
events through fragmentation of a 6 quark, as done in the standard PYTHIA, is too low 
(the ratio of the signal to the background is too small). For experimental feasibility studies 
of Be mesons, a very large sample of Be events is needed. Therefore BCVEGPY, which is 
powerful for generating B^ events only with full information in hadron-hadron collision, is 
very useful. 

Users may communicate with or give instructions to the program through an input file 
(input.dat). The output files include IsO.dat (the B^ total cross-section) or 3sl.dat (the 
B* total cross-section); pt.dat (the B^-pT distribution); shat.dat (the Bc-y/^ distribution); 
rap.dat (the Bc-Y rapidity distribution); pseta.dat (the Bc-t] pseudorapidity distribution) 
and grade.dat (the sampling importance function obtained by VEGAS). The input.dat file 
allows users to setup the generation parameters and requests, while the output files collect 
the relevant event information. 
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TABLE VI: The parameter values in the sequential order in the input.dat file. 

PMBC PMB PMC FBC 

PTCUT ETACUT ECM IBCSTATE IGENERATE IVEGASOPEN 

NUMBER NITMX 

NQ2 NPDFU NEV 

ISHOWER MSTP(51) 

IDWTUP MSTU(lll) PARU(lll) 

ISUBONLY SUBENERGY IGRADE 



The sequential order and the format of the parameters in the file input.dat should be the 
same as in the Table VI. The parameters specified in the input file are: 

• PMBC, PMB, PMC=: masses of the meson, h quark and c quark respectively (in units 
GeV); 

• FBC=: decay constant for the Be meson (in units GeV); 

• PTCUT=: pt cut for the Bc{B*) meson (in units GeV, value can be freely selected by 

users); 

• ETACUT=: Y cut for the Bc{B*) meson (value can be freely selected by the users); 

• ECM=: total energy for the hadron coUision (in units GeV); 

• IBCSTATE=: state of the 5^(5*) meson: IBCSTATE=1 for B^ll^So] and IBCSTATE=2 
for B*[1'S,]; 

• IGENERATE=: whether to generate complete events. IGENERATE=0, when users wish 

the simulation to stop after the generation of the 'final state' containing the B^. meson, 6-jet 
and c-jet of the subprocess gg Bc + b + c; IGENERATE=1, when users wish that complete 
events including the Be production are to be generated. In the latter case, IDWTUP=1; 

• IVEGASOPEN=: whether switch on/off the VEGAS subroutine: IVEGAS0PEN=1 for 
using VEGAS; IVEGASOPEN=0 for not using VEGAS; 

• NUMBER=: total number of times for calling the integrand (VEGAS parameter, see 
VEGAS manual). The parameter is needed only when IVEGAS0PEN=1 in each iteration; 

• NITMX=: upper limit for the number of iterations (VEGAS parameter, see VEGAS 
manual). The parameter is needed only when IVEGAS0PEN=1; 

• NQ2=: choice of Q^, the type of the characteristic energy scale squared in the production 
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(in units GeV^). Here seven choices are available: i). NQ2=1: — s/A {s is the squared 
center-of-mass energy of the subprocess); ii). NQ2=2: Q"^ — s; iii). NQ2=3: = p^^^ + 
m\- iv). NQ2=4: = {^ p ^^ + ml^ + ^pl, + ml + ^pl, + mlf; v). NQ2=5: = 
(v^Ptb, + + \lvlb + ^1 + v/ptc + "^c)V9; vi). NQ2=6: = + ^2 ^j^g -^^ 
parton distribution functions and in the coupling to the 6-quark line, and = 4m^ for the 
as in the coupling to the c-quark line; vii). NQ2=7: — p^f^ + ml; 

• NPDFU=: choice of the coUision type of hadrons. The assignments can be found in 
PYTHIA manual, e.g. NPDFU=1 for p - p and NPDFU=2 for p - p; 

• NEV==: number of the events for the hadronic production h + h ^ Bc + ■ ■ • {h means a 

hadron) to be generated; 

• 1SH0WER=: whether to switch on/off the showers, including initial and final states, 
multiple interactions, hadronization; e.g. ISH0WER=1 for 'on' and ISHOWER==0 for 'off' 
(see PYTHIA manual); 

• MSTP(51)=: choice of the proton parton-distribution set; e.g. MSTP(51)=2 for CTEQ3M; 
MSTP(51)=7 for CTEQ5L; MSTP(51)=8 for CTEQ5M etc. (PYTHIA parameter, see 
PYTHIA manual); 

• IDWTUP=: master switch dictating how the event weights and the cross-sections should 
be interpreted (PYTHIA parameter, see PYTHIA manual); e.g. when IDWTUP^l, parton- 
level events have a weight at the input to PYTHIA. Events are then either accepted 
or rejected, so that fully generated events at the output have a common weight; when 
IDWTUP=3, parton-level events have a unit weight at the input to PYTHIA i.e. they are 
always accepted; 

• MSTU(111)=: order of in the evaluation in PYALPS (a PYTHIA routine for calculating 
as, see PYTHIA manual); e.g. MSTU(111)=1 for leading order (LO); MSTU(111)=2 for 
next leading order (NLO); 

• PARU(111)=: constant value of ag (see PYTHIA manual), which is used only when 
MSTU(111)=0; 

• ISUBONLY=: whether to keep the information only of the sub-process gg Bc{B*)+b+c; 
ISUBONLY=0 for the full hadronic production, i.e. the structure functions are connected; 
ISUB0NLY=1 for the subprocess only. 

• SUBENERGY=: the energy (in units GeV) of the sub-process gg B^{B*) + b + c. It is 
needed only when ISUB0NLY=1; 



35 



• IGRADE=: whether to use the grade generated by previous running VEGAS when 
IVEGASOPEN=0; IGRADE=1 means to use; IGRADE=0 means not to use. 

In the package, two subroutines PHASE-GEN and VEGAS are needed when integrating 
over the phase space: 

The subroutine PHASE-GEN(YY,ET,WT) 

The subroutine is included in the file sqamp.for and is called by another subroutine 
PHPOINT in genevnt.for. The purpose is to evaluate the allowed phase-space points for the 
sub-process gg Bc{B*) + b + c when the center-of-mass energy of the sub-process is fixed, 
and to return a non-zero weight for each allowed phase space point. The energy-momentum 
conservation is integrated out in order to reduce the dimension of the phase space integration 
by four. Since the subprocess is a three body final state, the nine-dimensional phase space 
integration of the process turns into a five-dimensional one with a proper Jacobi determinant, 
when PHASE-GEN has been applied. 

The variables in the routine are: 

• YY(5)=: a five-dimensional random number with a range from to 1 for each dimension, 
corresponding to the five independent integration variables for the phase-space; 

• ET=: center-of-mass energy for the subprocess (in units GeV); ET can be chosen freely 
when ISUB0NLY=1, otherwise when ISUBONLY=0, ET is determined by PHPOINT; 

• WT=: the returned weight for each generated phase-space point. 

The subroutine VEGAS (FXN,NDIM,NCALL,ITMX,NPRN) 

• FXN=: the integrand calling the function TOTFUN in genevnt.for; 

• NDIM=: number of integration dimensions for the generator; NDIM=5 when 
ISUB0NLY=1; NDIM=7 when ISUBONLY=0; 

• NCALL=: maximum total number of the times to call the integrand in each iteration 
set by the user; 

• ITMX=: maximum number of allowed iterations set by the user; 

• NPRN=: print out level; (see VEGAS manual) e.g. NPRN=2, when printing out only 
the cross section values and errors. 

Note that the units of all the output data is well explained in the programme. 
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TABLE VII: Generation parameters used in the sample generation. 



INITIAL PARAMETERS 

Be IN l^^o 

GENERATE EVNTS 30000000 FOR TEVATRON ENERGY(GEV) 0.20E + 04 

M_{Bc}=6.400 M_{B}= 4.900 M_{C}=1.500 f_{Bc}=0.4800 

Q2 TYPE= 1 ALPHAS ORDER= 2 

PARTON DISTRIBUTION FUNCTION: CTEQ 3M 

USING PYTHIA MODEL FOR IDWTUP= 3 

PTCUT =0.000 GeV 

NO RAPIDITY CUT 

USING VEGAS: NUMBER IN EACH ITERATION= 300000 ITERATION= 27 
END OF INITIALIZATION 



C. Generator checks 

The whole Fortran package is checked by examining the gauge invariance of the amplitude. 
The matrix element vanishes when the polarization vector of an initial gluon is substituted 
by the momentum vector of this gluon. 

For further checking, we have performed several test runs. By setting the parameter 
ISUB0NLY=1 in the input.dat file, we obtain the transverse momentum and rapidity 
Y distributions for the produced Be, and the total cross-section for the sub-process gg — >■ 
Bc{B*) + 6 + c. The results, which are shown in the previous section, coincide well with 
several groups' results ^. 

When running the programme, the initialization is shown as a screen snap-shot in Ta- 
ble |VTT| Some output data are shown by Figs lllll'^ 

IV. CONCLUSIONS 

A Be meson generator BCVEGPY for hadronic collisions, based on the dominant sub- 
process gg Bc{B*) + b + c, has been developed and well-tested. The generator has 
been interfaced with PYTHIA, which takes care of producing the full event and filling the 



37 



standard PYTHIA event common block. In view of the prospects for Be physics at Tevatron 
and at LHC, the generator offers a valuable platform for further experimental studies. 
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APPENDIX A: THE HELICITY FUNCTIONS FOR THE AMPLITUDE 

In this Appendix, we show an example how to calculate the functions E^j^k- 

For evaluating the inner product {p-q), we introduce the notations k±, k± for a light-like 
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momentum k^^ and use the Weyl representation for 7- matrices: 









, (i = 1,2,3) 


(Al) 


, f = 






^1 









with Pauh matrices 



1 

1 



-i 

1 



1 
-1 



(A2) 



and 



(A3) 



k± = ko±h, k± = k.^ + iky = \k±\e''^'' = ^k+k_e''^ 

where 1 is the 2x2 unit matrix and cr* (i=l,2,3) are the Pauh matrices. Note that here 
we always have k+ > 0, A;_ > 0, because k+k^ = k^ — k^ — k^ + ky > due to k'^ — 
^0 ~ ~ — 0- By suitable choice of the phase we introduce the Weyl spinors 

/ \ / n \ 



\k 



k, 

"kZe'^^ 





, \k-) 






-^/kZ 



(A4) 



and 



{k^ ■ k2) = {ki-\k2+) = sJkZk^+e'^' - ^Jk,+k2-e'^' 



(A5) 



which appear to be explicitly antisymmetric. For the spinor product (A;i+||^3|A;2+), where 
A;f = 0(i = 1,2,3), 

{k,+ mk2+) - {k,+ \ks.){ks.\k2+) 

= — {ki+k2+h- - ki+k2±kl_^ - kl^k2+h± + k\^k2±h+) , (A6) 

and for the spinor product involving the polarization vector e{sz) of B*, 

qo+ 



(A7) 



where P' ^ P - ^^Qo- 



41 



Since the spinor products and the double inner spinor products are used frequently in 
the program, we take i/i {i — 1, ■ • ■ , 42) to denote all the possible spinor products appearing 
in the computation, e.g. yi = (/ci_|_|^^2ko+) and Xi,X2 to denote the double inner spinor 
products: 

Xi^ {qo- ki){qo- k2), X2 = {qo ■ h)* {qo ■ ki)* . (A8) 

To simplify the results, we have introduced the light-like momenta q^^, ((^^i Qci ?c2 ™ 
terms of time-like momenta qf,i, qb2, ?ci a-nd qc2 respectively: 

^qbi ■ Qo ■^1b2 ■ qo 

Qci = Qci - :r^qo , 4^2 = ?c2 - ;^^go ■ (A9) 
2?ci • go 25c2 • go 

When the functions Em,i,k ijn = l,---,5,8,9; k = l,---,64) are given, Em,j,k (j = 
2, ■ ■ ■ , 4) can be obtained by interchanging the initial gluon momenta and the final quark 
flavors. As shown below, due to the properties of the helicities, the functions may be simpli- 
fied, and may even become zero for some helicities. Here we list fi'm,!,! {tti — 1, • • • , 5, 8, 9) 
for an explicit example, 
4z/i 



^1,1,1 = —(jrib^iiyii - yi2)yi9y27 + yioiyaAVu - y^yl)) 

X\ 



+ yio(?/2o(?/2?/5 - yliyli) + T^cy^bylS) > 



—4 / 

^2,1,1 = (r«6^2/9 (2/112/192/27 + 2/IO2/34I/34) 

X\ 

+ yio{-(y9y2oy24y34) +^c^(y2y29y3i + 2/9^35^35)}) > 

-c/3,1,1 = ( (?«c 2/17 - yiiy32)2/i7 + {y5y5~y3m4) 

Xi 

+ y20{(y7 - 2/2)2/5 + 2/242/34} - ^c^2/352/35) > 

-£^4,1,1 = (2/9(2/ll2/322/?7 + ("^b^2/34 " 2/202/24)2/34) 

Xi 

+ rnc^{{y2 - 2/7)2/292/31 + 2/9(2/352/35 - yiryn)}) > 
_ 2yiyi6y27y29 

-C'5,1,1 — , 

Xi 

4 , 

-£^8,1,1 = — (y9(yi2 - yis) ("^6^2/34 - ^202/24)2/34 + 2gci • ^i(go • g6i)(^2 • /^i)* 2/272/362/35 + 

Xi 

"^c^{ (2/2 (2/12 - 2/18) - 2/72/16)^29^31 + y9((2/l2 - 2/18)^35^35 " 2/16^192/29)}) (AlO) 
£^9,1,1 = 0, (All) 

where the index c means charge conjugation of the spinor, \p_) — \p+y. 
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APPENDIX B: POLARIZATION VECTOR FOR B*[l^Si] MESON 



The expressions of the polarization vectors depend on the gauge choice. Here we choose 
a cartesian basis for the polarization vectors: 

1 



Pi - 



(P,,0,0,Po), 

P^Po 



P^M\JPi-Py'-P^ 



.^/Pi-Py'-P^, 



^Pi-Py'-Pf 



P'xPz 



^jPi-P'y-PV' 



( 



ey{P) 



PyPo 



.0, 



^JPS - P^ 



P P 



y-^ z 



(Bl) 



\^Pi-P^^/Po'-Py'-Pf ^Pi-P^-P'z ^Pi - P^^Pi - P'y - P! 

which satisfy the conditions 



(B2) 



e^-P-O, e^.e^ = -5*\ it,j=x,y,z) 



(B4) 



APPENDIX C: ROUTINES AND FUNCTIONS FOR THE HELICITY AMPLI- 
TUDE 

In this Appendix, subroutines and functions for calculating the helicity amplitudes for 
the sub-process gg Bc{B*) + 6 + c are explained. 
SUBROUTINE BUNDHELICITY 

Purpose: to compute the helicity amplitude of the bound state part, b + c ^ Bc{B*), where 
Be is the lowest state of I^Sq and B* is the lowest one of l^^*!, with the expressions for the 
polarization as presented in appendix B. 

Integer IBCSTATE-: state of the double heavy meson, IBCSTATE=1 for Bc[l^So]] 
IBCSTATE=2 for B*[l^Si]. 

Real*8 BUNDAMP(4)=: four hehcity amplitudes of the bound state part c + b^ Bc{B*)] 
BUNDAMP(4)=0 for c + 6 ^ B^; BUNDAMP(4)=1 for c + 6 ^ B^. 

Real*8 POLAR(4,3)=: four Lorentz components of the three polarization vectors ,{k — 
X, y, z) for the P* state. 
SUBROUTINE FREEHELICITY 
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Purpose: to compute the helicity amplitude of the process gg — > b + c + b + c. The functions 
^Fm ("^ = 1) ■ ■ ■ ) 5), Em,j,k = 1, ■ ■ ■ , 9; j = 1, ■ ■ ■ , 4; /c = 1, ■ ■ ■ , 64) are defined in the 
body of the paper. 

Real*8 PMOMUP(5,8)=: the momenta PM0MUP(5,j) (j = l,---,8) in the process: 
PM0MUP(5,1) for the gluon-1, PMOMUP(5,2) for the gluon-2, PMOMUP(5,3) for the 
Bc{B*), PMOMUP(5,4) for the 6, PMOMUP(5,5) for the c, PMOMUP(5,6) for the 6, 
PMOMUP(5,7) for the c, PMOMUP(5,8) for the go (the hght-hke reference momentum). 
Real*8 PMOMZERO(5,8)=: eight hght-hke momenta obtained from PMOMUP(5,8) ac- 
cording to Eq.(j2ni) respectively. 

Real*8 COLMAT(5,64)=: M'j^'^''^''^^'^"^'\Xi = ±, m = 1, ■ ■ ■ , 5). 

Integer IDP, IDQO, IDBl, IDB2, IDCl, IDC2, IDKl, IDK2=: symbols (codes) for the parti- 
cles in the processes: IDP for Bc{B*), IDQO for the reference massless fermion, IDBl for 
6-quark, IDB2 for 6-quark, IDCl for c-quark, IDC2 for c-quark, IDKl for gluon-1, IDK2 for 
gluon-2. 

Complex*16 YUP(42)=: 42 possible non-zero spinor products {yi, i = 1, 2, ■ ■ ■ , 42). 
Complex*16 XUP(2)=: two possible double-spinor inner products (xi,X2). 
Real*8 PROPUP(14,4)=: 14 possible denominators from the propagators appearing in the 
four basic groups of the Feynman diagrams {cb, be, cc, bb). 

Complex*16 BFUN(9,4,64)=: Em,j,k (m = 1, 2, ■ ■ ■ , 9; j = 1, ■ ■ ■ , 4; A; = 1,2, ■■■,64). 
SUBROUTINE AMP2UP 

Purpose: to compute the square of the helicity amplitude of the process gg Bc{B*) + b + c 
by connecting the amplitude of the bound state part to that of the free quark part. The 
helicities of the intermediate quarks c, b are summed up, the whole amplitude is squared, 
and then all the 16 possible squared helicity amplitudes (if particle polarizations in the final 
state are not measured) are summed up. Here we also consider the three kinds of color flows 
explicitly. 

Real*8 AMP2CF(3) =: three results of the square of the whole amplitude corresponding to 

the three kinds of color flows for the sub-process gg Bc{B*) + b + c. 

Real*8 FINCOL(5,16) =: 16 possible helicity amplitudes for the process, where the helicities 

of the intermediate quarks have been summed up. 

SUBROUTINE FIRST 
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Purpose: to compute the nine basic basic functions Em,i,k{jn = 1, ■ ■ ■ , 9; = 1, 2, ■ • • , 64). 
The correspondence of the nine functions to the Feynman diagrams is shown in Table |n] 
and TableEm In the present case: IDK1=1, IDK2=2, IDP=3, IDQ0=8, IDB1=4, IDB2=6, 
IDC1=7, IDC2=5. 
SUBROUTINE SECOND 

Purpose: to compute the nine basic functions, Em,2,k{jn = l,---,9;k = l,---,64) ob- 
tained from E^^i^k by gluon exchange. Here IDK1=2, IDK2=1, IDP=3, IDQ0=8, IDB1=4, 
IDB2=6, IDC1=7, IDC2=5. 
SUBROUTINE THIRD 

Purpose: to compute the nine basic basic functions, Efn^k{m = l,---,9;k = l,---,64) 
obtained from Em,i,k by 6, c quark exchange and 6, c antiquark exchange. Here IDK1=1, 
IDK2=2, IDP=3, IDQ0=8, IDB1=7, IDB2=5, IDC1=4, IDC2=6. 
SUBROUTINE FOURTH 

Purpose: to compute the nine basic basic functions, i?m,4,fc(^ = !;■■■) 9; fc = l,---, 64) ob- 
tained from E„,3,fc by gluon exchange. Here IDK1=2, IDK2=1, IDP=3, IDQ0=8, IDB1=7, 
IDB2=5, IDC1=4, IDC2=6. 
FUNCTION DOTUP(I,J) 

Purpose: to compute the dot-product of two momenta. 
Integer I,J=: codes of the two particles. 
FUNCTION INPUP(IP,JP) 

Purpose: to compute the inner product of two spinors with light-like momenta based on the 
formula presented in APPENDIX A. 
Integer IP, JP=: codes of the two particles. 
FUNCTION SPPUP(IP,KP,JP) 

Purpose(s): to calculate the spinor product with the formula presented in Appendix A. 
Integer IP,JP,KP=: codes of the three particles. 

FUNCTION POLSPPUP(I) 

Purpose: to compute the spinor product involving the B* meson polarization vectors with 

the formula presented in Appendix B. 

Integer I=: one of the three polarization vectors for B*. 
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FIG. 9: Bc-pT and Bc-Y (rapidity) distributions for the subprocess gg ^ Be + b + c with different 
center-of-mass energies \/l =20 GeV(I), 30 GeV(II) and 60 GeV(III). The sohd hne shows the 
present results and the dotted hne shows the resuhs obtained in Ref. 3|- 
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FIG. 10: Integrated cross sections for the subprocess gg Bc{B*) + 6 + c, with Ug = 0.2, mt, = 
4.9 GeV and rric = 1.5 GeV. The sohd (dotted) hne corresponds to the Be (B*) production. 
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FIG. 11: Bc-pT and \/l distributions for the CTEQ3M parton distribution function by using in 

the ncxt-to-lcading order (NLO) and adopting the characteristic energy scale squared = s/4, 
where s is the squared center-of-mass energy of the subprocess. 




FIG. 12: Bc-Y (rapidity) and B^-Yp (pseudorapidity) distributions for the CTEQ3M parton dis- 
tribution function, by using ag in the next-to-leading order (NLO) and adopting the characteristic 
energy scale squared = s/4, where s is the center-of-mass energy of the subprocess. 
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